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Abstract. In this paper we analyse the electronic properties of Dirac electrons in 
finite-size ribbons and in circular and hexagonal quantum dots. We show that due 
to the formation of sub-bands in the ribbons it is possible to spatially localise some 
of the electronic modes using a p — n ~ p junction. We also show that scattering, 
by an infinitely-massive wall, of confined Dirac electrons in a narrow channel, induces 
mode mixing, giving a qualitative reason for the fact that an analytical solution to 
the spectrum of Dirac electrons confined in a square box has not been found, yet. A 
first attempt to the solution of the square billiard is presented. We find that only 
the trivial case k = has a solution that does not require the existence of evanescent 
modes. We also study the spectrum of quantum dots of graphene in a perpendicular 
magnetic field. This problem is studied in the Dirac approximation, and its solution 
requires a numerical method whose details are given. The formation of Landau levels 
in the dot is discussed. The inclusion of the Coulomb interaction among the electrons is 
considered at the self-consistent Hartree level, taking into account the interaction with 
an image-charge density necessary to keep the back-gate electrode at zero-potential. 
The effect of a radial confining potential is discussed. The density of states of circular 
and hexagonal quantum dots, described by the full tight-binding model, is studied 
using the Lanczos algorithm. This is necessary to access the detailed shape of the 
density of states close to the Dirac point when one studies large systems. Our study 
reveals that zero energy edge states are also present in graphene quantum dots. Our 
results are relevant for experimental research in graphene nanostructures. The style 
of writing is pedagogical, hoping that new-comers to the subject can find this paper a 
good starting point for their research. 



PACS numbers: 73.21.Hb, 73.21.La, 73.23.-b, 73.43.-f, 73.63.Kv, 73.63.Nm 
1. Introduction 

Graphene was discovered in 2004 at the Centre for Mesoscopic and Nanotechnology 
of the University of Manchester, U.K., directed by A. K. Geim [T], [2j. Previously, 
graphene was known only as an intrinsic part of three-dimensional systems: as individual 
atomic planes within graphite or its intercalated compounds and as the top few layers 
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in epitaxially grown films [3]. In certain cases, it was possible to even grow graphene 
monolayers on top of metallic substrates and silicon carbide |3j. However, coupling with 
the substrate did now allow studies of electronic, optical, mechanical, thermal and other 
properties of graphene, which all became possible after individual graphene layers were 
isolated. There are by now a number of review papers on graphene available in the 
literature, both of qualitative [H [5], [H [3, [8] and quantitative [9l [10] nature. 

The original method of graphene isolation is based on micro-mechanical cleavage 
of graphite surface - the so called scotch tape method. This method, however, has a low 
yield of graphene micro-crystallites. Recently, a new method [121 [H], based on liquid- 
phase exfoliation of graphite has proved to produce a large yield of graphene micro- 
crystallites, with large surface areas. A chemical approach to graphene production has 
also been achieved using exfoliation-reintercalation-expansion of graphite [T3] . 

It is by now well known that graphene is a one-atom thick sheet of carbon atoms, 
arranged in a honeycomb (hexagonal) lattice, having therefore two carbon atoms per unit 
cell. The material can be considered the ultimate thin film. In a way, this material was 
the missing allotrope of pure carbon materials, after the discovery of graphite [HI [TS] . 
diamond [IE], fullerenes, and carbon nanotubes [H]. In fact we can think of graphene 
as being the raw material from which all other allotrope's of carbon can be made [HIH]. 

Although the Manchester team produced two-dimensional micro-crystallites of 
other materials [2], graphene attracted a wide attention [I7] from the community due 
to its unexpected properties, both associated with fundamental and applied research. 

Graphene has a number of fascinating properties. The stiffness of graphene has 
been proved to be so large, having a Young modulus E ~ 1.0 TPa, that makes it 
the strongest material-stiffness ever measured [HI [18]. In addition, the material has 
high thermal conductivity [T9], is chemically stable and almost impermeable to gases 
[20] . can withstand large current densities [1], has ballistic transport over sub-micron 
scales with very high mobilities in its suspended form (/i ~ 200.000 cm^-V~^-s~^) 
[H [211 [22], and shows am-bipolar behaviour [1]. Ballistic transport is in general 
associated with the observation of conductance quantisation in narrow channels [23) . 
which have recently been observed [24]. The above properties and its two-dimensional 
nature makes graphene a promising candidate for nano-electronic applications. 

From the point of view of physical characterisation, we are interested in the 
mechanical properties of graphene, its electronic spectrum, its transport properties of 
heat, charge and spin, and its optical properties. The high stiffness of the material 
is responsible for the micro-crystallites to keep their planar form over time, without 
rolling up, even when graphene is held fixed by just one of its ends [H]. Interestingly, 
graphene is now being used, by combining electrostatic deposition methods and the 
chemical nature of the surrounding atmosphere, to produce rolled up nanotubes with 
controlled dimensions and chiralities [25]. The electronic spectrum of graphene and 
of graphene bilayer has been measured by angle resolved photoemission spectroscopy 
(ARPES) [26l [27l [2HI [29] , fitting well with tight-binding calculations using a first nearest 
neighbour hopping t ~ 3 eV and a second nearest neighbour hopping t' ~ 0.13 eV [30]. 
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The transport of heat has been measured experimentally and studied theoretically 
in few papers |T9l [HI [32l [HH [34|, and more research is needed to fully understand 
its properties, especially since the irradiation of graphene with laser light has been 
found to heat up locally the system [19]. The transport properties of graphene on 
top of silicon oxide have been extensively studied, but in a suspended geometry the 
few available experimental and theoretical studies are still recent [22l [35l [36l [37]. Of 
particular interest is the contribution of phonons to the transport properties of graphene. 
Although in the beginning of graphene research phonons were considered not important, 
recent experimental results show that this is, in fact, not the case [211 [38]. 

Of particular interest is the finite conductivity of graphene at the Dirac point, 
aD, whose measured value is in contradiction to the naive single particle theory, which 
predicts either infinite or zero resistance. The value of ct/) is of the order of 

CTD ^ A4^ , (1) 

with A a number unity order. We emphasise that Eq. ([T]) is the value for the conductivity 
of the material at the Dirac point, and not the conductance of a narrow channel. 
On the other hand, there is a discrepancy between the more elaborated theoretical 
descriptions of an and the experimental measured values, since the theory predicts the 
value a^"^"^' = a/j/vr [39l [40l [HI [42] , and most of the experiments measure a value given 
by Eq. ([T]). Adding to the problem, two experimental groups reported measurements 
of the conductivity of graphene consistent with the theoretical calculation [43|,[44]. In 
this context, it should be stressed that whereas the result of ajj obtained in Ref. [39] 
comes about due to an increase of the density of states due to disorder (albeit small) at 
the Dirac point, the value for an computed in Refs. [40l [41] is based on the existence 
of evanescent waves in clean graphene ribbons with large aspect ratio W/L (W is the 
width and L in the length of the ribbon). The recent experiments [43l [44] seem to 
confirm this latter view of the problem, since the value cr^*^"^- is only measured in the 
regime W/L S> 1. The transport of spin in graphene was studied experimentally in few 
publications [45|, [46|, [471 [M], and much work remains to be done. 

The optical properties of graphene and of bilayer graphene have only recently been 
studied experimentally [49l [50], in contrasts to the corresponding theoretical studies. 
The first theoretical study of the graphene's optical absorption was done by Peres et 
al. [39], followed by several studies by Gusynin et al. and reviewed in [5l]. The most 
relevant aspect was that the infrared conductivity of graphene, for photon energies larger 
than twice the chemical potential, has a universal value given by [39l [5ll [52] 

Also for the bilayer, the theoretical studies preceded the experimental measurements 
[55l [56l [57] . Due to the existence of four energy bands in bilayer graphene, its optical 
spectrum has more structure than the corresponding single layer one. 

It was experimentally found that for photon energies in the visible range [53] Eq. 
([2|) also holds within less than 10% difference [54]. This result makes graphene the 
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first conductor with light transmissivity, in the frequency range from infrared to the 
ultraviolet, as high as 

T ~ 1 - 7ra ~ 98% , (3) 

with a the fine structure constant, in a frequency range that extends from the infrared to 
the ultraviolet. This makes obvious that graphene can be used as a transparent metallic 
electrode, having found applications in solar cell prototypes [581 El] and in gateable 
displays [60]. The same conclusions are obtained from studying the optical conductivity 
of graphite [H]. The transparency of graphene has an obvious advantage over the more 
traditional materials, used in the solar cell industry, indium tin oxide (ITO) and fluorine 
tin oxide (FTO), which have a very low transmission of light for wave-lengths smaller 
than 1500 nm; in the visible range the transparency of these two materials is larger 
than ~ 85%. Furthermore, these traditional materials have a set of additional problems 
[Ml [60], such as chemical instability, which are not shared by graphene. On the other 
hand, ITO and FTO have a low resistivity (~ 5f2-m), a figure that graphene cannot 
match, if one leaves aside the possibility of graphene films deposit from solution. 

Finally, the interaction of graphene with single molecules allows to use graphene 
as a detector of faint amounts of molecules [62j and to enhance, in a dramatic way, 
the sensibility of ordinary transmission electron microscopes [I2l [63|, allowing the 
observation of adsorbates, such as atomic hydrogen and oxygen, which can be seen 
as if they were suspended in free space. 

Many of the above properties are expected to be present in graphene nanoribbons. 
On the other hand, aspects related to the quantum confinement of electrons in graphene, 
both considering confinement in one- (quantum wire) or two- (quantum dot) dimensions, 
is expected to bring new interesting phenomena. The present state of the art of 
material manipulation technologies does not allow to produce structures, on a top- 
down approach, smaller than 10 nm. Nevertheless, many aspects associated with the 
quantum confinement of electrons in graphene either due to narrow constrictions or 
due to the formation of quantum dots have already been investigated experimentally 
[Ml iSl Ei, |67l iH and theoretically [69l [70l [HI [72]. One important consequence of 
quantum confinement is the appearance of an energy gap in the electronic spectrum of 
graphene, an important characteristic if graphene is to be used as a material to build 
nano-transistors [73]. In this paper we will address several properties of confined Dirac 
electrons, by considering both nano-wires and quantum dots of graphene. In doing 
this we use both the continuous Dirac approximation and the tight-binding description, 
choosing which one is more appropriate to the given problem. 

2. The tight-binding model and the Dirac approximation 

The experimental work cited in the introduction constitutes a vast evidence that the 
low-energy theory of electrons in graphene is described by the two-dimensional Dirac 
equation, which is obtained as an fc-p expansion around the Dirac points in momentum 
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Figure 1. (colour on-line) An honeycomb ribbon, with zigzag edges (top and bottom) 
and arm-chair edges (vertical ones), with the carbon atoms belonging to the sub- 
lattices A and B clearly differentiated from each other (with several carbon atoms 
represented). The lattice unit vectors Oi and a2 are also shown. 



space [9l[74]. In Figure [H we represent a finite-size ribbon of an hexagonal lattice. Two 
features are of importance: the first is that the lattice is not a Bravais lattice, being 
instead made of two inter-penetrating triangular lattices, giving rise to two geometrically 
nonequivalent carbon atoms, termed A and B; the second is that there are two different 
types of edges present - zigzag and arm-chair edges. These types of edges play a different 
role in the physics of the ribbon. In particular, note that the zigzag edges are constituted 
by a single type of atom, B on top and A at the bottom (in the case of this figure). As 
is shown in Fig. [H we can choose the direct lattice vectors to be the following: 

ai = |(-l,v^), (4) 

a2 = ^(l,v^), (5) 

where a = 2.46 A is the lattice vector length. As a consequence, the reciprocal lattice 
vectors are 

27r , 1 , , . 

27r, 1 , 

The so called Dirac points in the honeycomb Brillouin zone are conveniently chosen to 
be 

/f=5j(1.0). (8) 

/f'=|(-1.0). (9) 

If one considers only the hopping process t (first nearest neighbour hopping), the tight- 
binding Hamiltonian is very easily written as (Nc is the number of unit cells in the 
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solid) 

i/=-t^(4A- + H.c.), (10) 

i,cr 

where aj^ creates an electron with spin projection a in the vr-orbital of the carbon atom 
of the sub-lattice A, and of the unit cell i; a similar definition holds for bl^. The exact 
diagonalisation of this problem is straightforward, leading to 



E± = 3 + 2 cos(aA;a;) + 4 cos 2 """^^ ^ 2^ ^^ 

We can expand this relation near the Dirac points, obtaining (fc = q + K) 

E± ~ ±vf\c^\ , (11) 

which is a massless Dirac-like linear dispersion relation, where the velocity of light 
is substituted by vp = ^^t ~ 10^ m/s, the Fermi velocity. To obtain the effective 
Hamiltonian obeyed by the electrons near the Dirac points we write the matrix 
Hamiltonian in momentum space as 

with Sfc given by (fc = g + K) 

Sk = l + e'^'^'e''^''' + e'^-'^^e'"-''^ 

^{(lx + ^qy), 

leading to the effective Hamiltonian (one valid near K and the other near K' ) 

HKiq) =VF(T*-q (13) 
HK'iq) =Vf(t ■ q, (14) 

with cr = {ax,(Jy) and cr* = (cr^,— cTj,). The Hamiltonian (fT3|) . or alternatively (fT4l) . 
will be the starting point of our discussion. From Eqs. ( fT3l) and (fT4l l we can write the 
second-quantised Hamiltonian for electrons in graphene 

H ~ -ivF j dx(ly{i!\{r)(T ■ V^i(r) + ^5(r)cr* ■ V^2(t')), (15) 

where = (a|, 6j) (for i = 1, 2). 

Let us assume that it is possible to create a potential such that a term of the form 

^ = XI ^'F^(^ia(^i,'^ - blak<r), (16) 

is added to the Hamiltonian (fTOl) . In terms of the formalism used to write the effective 
Hamiltonian (fT3l) and ffT4l). this term is rewritten as 



V = vlma, , (17) 
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and it corresponds to the presence of a mass term in the Hamiltonian. This type of 
term can be generated by covering the surface of graphene with gases molecules |75j or 
by depositing graphene on top of boron nitride [76l [77l [78]. The eigenvalues of H^^ + V 
are easily obtained, leading to E = Vpff\k\'^ + m^Vp, with \k\ = ^Jk'^ + ky, and the 
same for Hk' + V . 

3. Confinement of Dirac electrons on a strip 

Our goal in this section is on deriving a mathematical framework for describing the 
effect of confinement on Dirac electrons. The confinement can be produced either by 
etching, by the reduced dimensions of the graphene crystallites, or by the application 
of gate potentials (here the Klein tunneling poses strong limitations on the use of such 
method). 

3.1. Boundary conditions and transverse momentum quantisation 

The mathematical description of the confinement requires to impose appropriate 
boundary conditions to the Dirac fermions. Although, for graphene ribbons, the two 
types of edges discussed above impose two different types of boundary conditions [82] . 
we shall use here the infinite mass confinement proposed by Berry and Mondragon 
[83] . For large ribbons, there will be no important difference between the two types of 
boundary conditions [84], except that the infinite mass boundary condition is not able 
to produce edge states [85], which are present in ribbons with zigzag edges. 

We shall generalise some of the results of Refs. [40l [83] by considering the case of 
Dirac fermions with a finite mass. The mass profile in the transverse direction [y) of 
the strip is represented in Fig. [3 The boundary conditions the wave function has to 
obey, at the spatial point where the mass changes from m to M, are derived considering 
the reflection of the wave function at the boundary. Let us first consider the refiection 

X It is important to comment here on this particular choice for the boundary condition. Using the 
k ■ p approach of DiVincenzo and Mele [74] one learns that vp 1/m, where m is the bare electron 
mass. On the other hand, one considers that graphene electrons can not propagate in a region where 
the material is absent, and therefore have zero velocity there. Due to the proportionality vp oc 1/m, 
this can be achieved taking the limit m ^ oo. Therefore we could think that confinement of Dirac 
fermions could be achieved with a position dependent Fermi velocity vpiu), that goes to zero at the edge 
of the strip. Unfortunately this program does not work. Berry and Mondragon boundary boundary 
condition[83j corresponds to a change in the nature of the spectrum, and is somewhat artificial in 
what concerns graphene. The consequences of the different boundary conditions are: the properties 
of the wave function at graphene edges do depend on the different choices of boundary conditions, 
but the bulk behavior of the electronic states is essentially the same for all them; very close to the 
neutrality point the choice of boundary condition do again matter, but at finite doping this is not the 
case anymore. The choice of Berry and Mondragon boundary condition[83] does introduce a certain 
degree of simplicity in the calculations. 
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Figure 2. Scheme of the mass confinement (along y) with mass m inside the strip 
and mass M outside. 



aX y — L/2. The wave function in the central region (7) is given by 



with Ok = aj:ci&n.{ky / kx) 



E — mVp 



2nA 



In zone II {y > L/2), the solution has the form 

i^ii{x,y) = T 



1 

fii(E) 



with 



and 



fn{E) 



Qy = ±1 



E-Mvl 

VfK^x - iqy) ' 



lE-_-Mh^ _ 



^ik^x^ (18) 



(19) 



(20) 



(21) 



(22) 



where the energy values are given by the same expression as that for zone /. As we 
want to take the limit M — > oo, we will assume Mv\ > E, which implies that 



(23) 
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and thus 



fniE) 



E-Mvl 



E-Mvl 



(24) 



where the sign ± in front of the square root applies to the wave function that is 
propagating in the positive/negative y direction.. Imposing the boundary condition 
associated to the Dirac equation for a reflection at y = L/2, 



one obtains 



ih 



1 



V'ji 

Taking the limit M ^ oo the boundary conditions reduce to 



-L/2 



y=L/2 



+ 1, 



1 . 



(25) 



(26) 



(27) 



(28) 



Now one wants to write down the wave function of electrons propagating on the 
strip taking into account the confinement due to the mass term. The most general wave 
function is of the sum of two counter-propagating waves in the y direction 



where 



i^ix,y) 



x{y) = A 



Jkyy 



fiiE)e^'' 

It is always possible to redefine B such that 



+ B 



fi{E)e 



-ikyp 



-iSh 



(29) 



(30) 




x{y)=A\ + 5 ^ e-''\ (31) 



a procedure that proves useful later on. Imposing the boundary conditions ( l28l) . one 
obtains (considering the strip to be in the range < ?/ < L for simplicity) 

for the relation between the coefficients, and 

^ _ 1 - fi{E)e^''^ 1 + fi{E)-^e^'- 

1 - fj{E)-^e'o^ l + fi{E)e^'>' ^ ' 

for the energy quantisation. It is clear that the admissible values of ky are energy 
dependent. Considering the limit M — > oo one obtains, from Eq. (I32ll . the simpler 
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result A = B, and, from Eq. ( l33l ). the condition e* 
transverse momentum quantisation rule |83| 



71 mr 



where n = 0, ±1, ±2 



-1, which leads to the 



(34) 



Putting all together, the obtained results are summarised as: 



A 




( 1 ^ 




( Zn,k 1 












^ 1 ) 





Xn,k{y) = A 

where we have used k = k^, qn = ky^, s = ±1 = sign[£'] and 

k + iqn 



(35) 



Zn,k 



se 



If we ignore questions of convergence, we recognise that this form for 2;„ ^ does not 
require k or g„ to be real. We are only assuming k"^ + q^ > 0. The dependence of the 
energy E on k shows a number of sub-bands separated by energy gaps; this is shown in 
Fig. ^ for a ribbon 10 nm wide. It is clear that for such a narrow ribbon one has large 
energy gaps between two consecutive sub-bands. 

E(eV) 











-3 -2 -1 


rsj^ 1 2 3 


^^^^^^^^^^^^^^^^^^ 





kx (nm"') 



Figure 3. (colour on-line) Energy levels of a ribbon 10 nm wide. The transverse 
modes range from n = to rt = 8, for both particles and holes. 

Let US represent Eq. (l35l) as \'^n,k) = \Xn,k)G^^^, it is then simple to show that 
{^m,k'\'^n,k) = and that the normalisation coefficient A reads A = 1/{2^/L). Note 
that if on the strip we have a non-zero scalar potential IV, we will just have to substitute 
E hj E — V and replace k hj k, with k given by 



k' 



VTpTTl 



(36) 
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Figure 4. Representation of a step potential in a strip with lateral confining infinite 
mass. Zone I is for 2: < and zone II is for a; > 0. 



3.2. Dirac fermions in a strip with a step potential 

Let us now consider the scattering of Dirac fermions in a strip by a simple step potential, 
as represented in Fig. [H In the zone I one has V = 0, the wave function is given by Eq. 
(I35II : in zone //, with V > 0, the wave function is also given by Eq. (l35l) making the 
replacement k ^ k. If, in general, the step rises up at x = X the boundary condition 
takes the form 

i^n,k{X, y) + r„^„ _fc(X, y) = tni/j^-^iX, y) . (37) 
Solving for r and t gives 



r„ = 'll^^I^e'^^'' , (38) 

Equations (l38ll and (l39l l represent the reflection and the transmission amplitudes, 
respectively, for the transverse mode n. It is now a simple matter to compute the 
tunnelling transmission for an arbitrary conflguration of finite potential steps by using 
these two results combined with a transfer matrix method [78]. A particular case of this 
situation is the transmission through a potential barrier, for which is given in Fig. 



3.3. Trapped eigenmodes 

In this section we will show that it is possible to trap massless Dirac electrons in a ribbon 
of finite width L by creating a p — n — p junction [79]. The effect exploits the fact that 
the spectrum of a finite ribbon exhibits energy gaps. A similar study was done in Ref. 
[86] for the bulk case. In this case, such a confinement is possible for certain incident 
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angles of the eigenmodes on the potential walls [86]. The potential profile considered 
is shown in Fig. [H We show that the trapping of the eigenmodes requires evanescent 
modes in the x direction. This can be accomplished using a scalar potential. 

In order to solve this problem let as again consider the case of a potential step as 
in Fig. [6l In region I {x < w) the wave function has the form (1351 ) and in region // 
{x > w) the form would be the same with k replaced by 



One now makes the observation that if obeys the condition 

{E - V f 2 

vlK' ^ ^" ^ ^^^^ 
one has a propagating wave in region / and an evanescent wave in region II. Within 
the validity of Eq. (|4T]1 it is more transparent to write the wave function in region // 

as 



2^1 



IS 



IS' 



+ 



1 



with 



k 



la 




E-V 

vph 



(43) 



where s' = sign[i? ~ Of course, the same type of analysis holds if one had considered 
the step at X = 0, starting at the interface between regions /// and I. The trapping 
"mechanism" uses this fact. The wave function in region I of Fig. [6] is taken as a sum 
of two counter-propagating waves along the x direction, whereas in regions // and /// 



T(?„) 




Figure 5. Transmission coefficient, T{qn) — |tnP, through an energy barrier, of length 
w — 100 nm and height V — 50 eV, as function of the transverse quantisation quantum 
number n. The energy of the electron is taken as -E = 0.1 eV and the width of the 
ribbon is L — 500 nm. 
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Figure 6. Scheme of the confinement (along y) in a strip where a scalar potential 
well, of width w, was created. On the left one has an upper view, and on the right one 
has a side view. 

Table 1 . Values of w (in nm) for a given momentum g„ (in 1 /nm) . The parameters 
are s = 1, s' = -1, vp = 10^ m/s, E = 0.1 eV, V = 0.15 eV, and L = 500 nm. 

n w n w n w n w 

5 -59 6 -65 7 -78 8 -113 



only evanescent waves exists. Imposing the boundary conditions at a; = and x = w 
{w the width of the well) one obtains after a lengthy calculation the condition of the 
energy of the trapped eigenmodes 

sm{kw)F{E, V, Qn) + cos{kw)G{E, V, gn) = , (44) 

with 



~ 16 



and 



G{E, V, qn) — —{Zn^k + Ki,k + ^n,k + {^n,kf)i.^n,a 



(45) 



(46) 



Both F{E, V, qn) and G{E, V, qn) are pure imaginary numbers, as long as condition (1411 ) 
holds true. In order to give a flavor of the numerical solution of Eq. (I44l l we present 
its numerical solution in Table [H We have chosen the strategy of fixing the energy and 
looking for the values of w that satisfy Eq. (1441) . for different values of qn- 



4. Inducing mode mixing by scattering at a wall 

In this section we want to discuss the scattering of Dirac electrons when they propagate 
along a semi-infinite narrow channel and scatterer back at the wall located at the end of 
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the channel. This problem is intimately related to the possibility of finding a solution 
for the eigenmodes and eigenstates of trapped Dirac electrons in a square box. Our 
analysis hints at the reason why this solution has not been found yet. 



4.1. Definition of the problem 



Let us now consider massless Dirac fermions confined in a semi-infinite strip: a; < and 
< y < L. We will look for scattering states produced by the scattering at the wall due 
to an incoming wave from x — > —00. As before, the fermions are confined by an infinite 
mass term outside the strip. 

The scattering state is a sum of an incoming wave, with energy E, longitudinal 
momentum k, and transverse momentum g„, with a superposition of all the possible 
outgoing channels with refiection amplitude r„^^, and it can be written as 



^i{x,y) 



m=0 




(47) 



Since we are considering elastic scattering is an eigenstate), we must have 



I.e., 



We must distinguish two situations: 



(48) 



V e^^, (49) 
ql>E':kr^^iy^Z:^. (50) 

The second case corresponds to evanescent modes. We must choose this solution if we 
want the wavefunction to be convergent for x — > —00; for real km the choice of sign 
refiects the fact that we have only one incoming mode. We now simplify the notation, 
using the fact that k and n are fixed, and define 

k + iqn 



The parameter Zn is just a phase, \zn\' 
but for evanescent modes \zmf 7^ 1- For q^^ E"^, we obtain 



(51) 
(52) 

1; Zm is also a phase for propagating modes, 



E 
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4.2. Calculation of the reflection coefficients 

The boundary condition at the end of the semi-infinite strip a; = and < y < L, is 
(see Sec. [5] for details) 

^i(0,y) + z^2(0,y) = 0. (53) 

Note that the boundary condition at an infinite-mass vertical-wall is different from 
that of the horizontal case discussed before. Applying the boundary condition (l53l) to 
^n,k{.x,y) and after rather lengthy algebra (where the replacement — (m + 1) = m' is 
made at some stage and m' is redefined as m afterwards) one arrives at the condition 

00 

I{n,m)e'''^y = 0, (54) 

m=— 00 

with J(n, m) given by 

J(n, m) = [(1 + iZjn) ^n,m + ?^n,m(l + iZm)] 9(771 + 1/2) 
+ [(^ 

+ i)]e{-l/2-m) . (55) 

The crucial step in the derivation is the observation that the set of function's {(pm = 
g«9mi/^^ = 0, ±1,±2, ...} is overcomplete. In fact, the set of states with m even (or 
with m odd) is, by itself, a complete orthogonal set for functions defined in the interval 
< y < L. If follows that we obtain an equivalent set of conditions to Eqs. (I54ll by 
taking inner products of Eq. (l54ll with 0^ with m even (or m odd). We recall the inner 
products (choosing p even) to be 



1 /-^ . I ^' 

- rfye"*^'^''"^'"^^ = < 0, mis even, m 7^ p . (56) 

^ Jo 2 _ • „ J 1 



mis odd 



The fact that the integral (l56l l is not a Kronecker symbol shows that, in this case, the 
basis is overcomplete. Using Eq. (l56ll . we obtain 



2 

J(n,p) + V- —l(n,m)=0 p = 0, ±2 ± 4, . . . (57) 

^-^ i[v — m)T[ 

rriodd 

which is the central result of this section. This is the set of equations needed to calculate 
the r„ m coefficients. Naturally the convergence of the sum in ([57l) critically depends on 
the behaviour of rn,m with m. 

We could also have formulated the scattering problem a bit more generally, 
considering, for example, the case of a wave that approaches a wall at x = —D coming 
from X = +00. Naturally the modifications relatively to the solution found before 
cannot be much, given the symmetry of the problem. The first thing to note is that the 
boundary condition is slightly changed, being given by (see Sec. [5] for details) 

^,{-D,y)-t^2{-D,y) = 0. (58) 
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Working out the problem along the same lines as before, one learns that the final result 
can be obtained from the previous solution upon the replacements 

r r piD{k+km) /irq\ 

' n,m ' ' n,m^ ; \'~'^ ) 

Zm ^-l/Zm, (60) 

Zm ^-l/Zm, (61) 

where the transformation ( l59ll is obtained using the generator of translations, T(xo) = 
^ixop/h p ^\^Q momentum operator), and follows from the new position of the wall. 

The transformations (l60l) and fl6Tl) follow from the difference in the boundary condition 
between a right and a left vertical wall. One should note that the transformation (l59D 
is not a phase for evanescent waves. Using transformations (l59ll -(l6Tll in Eq. (l55l) one 
obtains 



I{n, m) [{1 + i/zm) 6n,m + e^^('^+'='"V„,„(l + i/zm)] e{m + 1/2) 

- [{1/z^m-l + Sn,-n.-l + e*^(''+'^— V„,_™_i ( l/z_^_l + z)] 0(-l/2 - (fi 

Naturally, the exact solution of the scattering problem rests upon the possibility 
of solving exactly the set of linear equations (l57ll . This task seems out of reach at the 
moment. The second approach is to solve numerically this set of equations. This leads 
to the conclusion that the summation over m has to be truncated at some value. To 
be concrete, let us consider the particular case of an incoming mode with transverse 
quantum number n = 0, such that the value of the incoming longitudinal momentum 
k originates Np propagating modes above the mode n = (0 < A'p < A^). Then the 
numerical solution of the problem has to satisfy the conservation of the probability 
density current (it must be one in this case) and the retained coefficients after the 
truncation have to converge upon increasing A^. As we show below, both these two 
conditions are satisfied for small A^. For the numerical solution we choose n = 0, 
go = 7r/2, such that the value oi k = 27i leads to n = and n = 1 as the only two 
propagating modes. We use a set of units such that L = 1 and hvp = 1. Further we 
take E > which implies that s = 1, since the scattering being elastic can not excite 
hole states, which have E < 0. The equations to be solved numerically are given in 



Appendix A 



The probability density current transported by the mode n is defined as 

Snix, y) = VF^jiix, y)(T^ljn{x, y) . (63) 

Since the motion is transversely confined, what is needed is the probability density 
current along the x-direction, which reads 



Sn{x)=VF dyi)l{x,y)a^'4)n{x,y) . (64) 

JO 

Using definition f l64ll . the total reflected flux density, S'^''^'^^\ obeys the sum rule 

Nv n 

QT,refl i |2 , ST^ '^'^^ Prn \ ,2 ^ 

= rO,0 ^rO,»n =1, (65) 

^ COS^O 
m=l 
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with 



arctan — , 



(66) 



with both Qn and kn real. In agreement with our expectations, the sum-rule is better 
fulfilled the larger N is, although modest values of do a good job as well. In fact, in 
the right panel of Fig. [7] one can see that the sum rule is fulfilled even considering only 
one evanescent mode (A^ = 3). 




0.5 
k/(2K ) 



Figure 7. (colour on-line) Left panel: Energy levels for n = 0, 1, 2, and k — 27r. This 
leads to Np = 1. Right panel: sum rule l|65p for Np = 1 as function of iV = 3, 4, . . . , 10. 
We have depicted only odd values of N, such that the total number of odd and even 
terms is the same (n = is considered even) , but our results are independent of this 
choice. 



In Fig. [8] we study, for the particular mode occupation defined in the left panel of 
Fig. [7l the evolution of the coefficients ro,„ as function of A^. We write each coefficient 
rQ p as rQ p = | ro,p | e*"'^-*' . In Fig. [H we plot the square of the modulus of rg^p (left 
panels) and the corresponding phase ao,p (right panels), separating the cases for which 
p = 0, 1 (propagating modes), which are represented in the two top panels, from those 
where p > 2, . . . 7 (evanescent modes), which we represent in the two bottom panels. 
A beautiful result emerges from this study. The wall introduces mode mixing and 
generates evanescent waves, whose contribution to the total wave function diminuishes 
upon increasing p. This result is quite different from that for Schrodinger electrons, 
where no mode mixing takes place. The fundamental reason is due to the fact that 
for Schrodinger electrons the transverse wave-function is the same for the incoming 
and outgoing waves. For Dirac electrons, on the contrary, the spinor of the incoming 
and outgoing waves change due to its dependence on either the incoming or outgoing 
momentum. 

Had we tried to force the solution of the problem using only one incoming and 
one outgoing propagating modes, with the same k value, and we would have obtained 
the trivial solution k = 0. This statement is easily proved as follows: we make the 
assumption that the total wave function should be a sum of two terms of the form 
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Figure 8. (colour on-line) Top panels: evolution of |ro,oP and |roap as function 
of N (left); evolution of the phases q;o,o and q;o,i as function of N (right). Bottom 
panels: the same as before but now for |ro,pp (left) and for ao.p (right), considering 
p = 2,...,7. 




+ 



+ 



1 



(67) 



where the phase-shift 5 was introduced. Let us now impose the boundary condition (1531) 
on the wave function (l67|l . Working out the calculation, one obtains two conditions that 
must be fulfilled simultaneously 

cos 5 ± sin(/5„^fc — 5) = , (68) 

with (3n,k defined from Zn^k = e*'^" * . It is clear that the two conditions in Eq. (l68l) cannot 
be satisfied, in general, at the same time, which precludes the proposal of Eq. f l67l) as a 
solution to the problem. In fact, the two conditions given by Eq. (l68l) are equivalent to 



5 



TT 



+ fTT A 5 = + £7r , £ = 0,1,2,..., 

TT, a situation that occurs only if 



(69) 



(70) 



which can only be true if 2f3n,k 

In 71" 

arctan — = — , 
k 2 

which finally is true only in the trivial case /c = 0. This means that the wave function has 
no X dependence and that the electronic density is p{x,y) = "^l^ f,{x,y)'^n,k{x,y) = l/L, 
constant everywhere (for finite m the density does show oscillations inside the box [80]). 
The exact solution of the square billiard with infinite-mass confinement is therefore a 
quite elusive problem |83| . 
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5. Confinement of Dirac fermions in quantum dots 

Let us now consider the confinement of Dirac fermions in quantum dots. Naively one 
would expect that the rectangular dot would have a simple solution (as it has in the 
Schrodinger case), since the wave function of the confined Dirac electrons (by an infinite 
mass term) in a strip can be written in terms of elementary trigonometric functions. 
In fact this is not the case. The only known case so far of an integrable Dirac dot 
(billiard), subjected to the infinite mass confinement, is the circular one. In order to 
solve this problem one needs the boundary condition obeyed by the wave function at 
the dot boundary. This was worked out by Berry and Mondragon [83] and the geometry 
they used is represented in Fig. [9l They considered a quantum dot represented by a 
domain D of arbitrary shape, separated by an outside region that we denominate OD. 
The boundary of the domain D is parametrised by an length arc s(a), where the vector 
normal to the surface of the dot at s is given by 

n(s) = cosq;(s) + sina(s) Cy (71) 

Imposing the condition of zero fiux perpendicular to the wall of the dot one obtains 




Figure 9. (colour on-line) Domain D with the boundary parametrised by s. Figure 
adapted from Ref. ^83j. The incident and reflected wave at s are both shown. 



(72) 



The constant B is determined working out the study of a refiecting wave at the boundary 
of the dot, when the mass in the region OD obeys the condition M ^ oo. The final 
result is i? = 1 [83], and the detailed calculation can be found in [Appendix B 



5.1. The circular dot with zero magnetic field 

Let us first write the free solutions of the Dirac equation in polar coordinates r and ip. 
In this coordinates the Dirac Hamiltonian and the wave function read [871 [88ll89ll90l [9T] 
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and 

respectively, where Jm{x) is the Bessel function of integer order m. A detailed derivation 



of these results is given in |Appendix C[ In order to obtain the eigenvalues of the electrons 
in the dot one has to apply the boundary condition (rf2l) . Note that because one has a 
circular dot a{s) = ip. This latter property makes it possible to satisfy the boundary 
condition (l72ll with the wave function (Tfil) alone. In fact, for a dot of radius R, one has 

siJ„+i(i?A;)e*("^+^)^ = JURk)e'"'^ie''' ^ sJ^+i{Rk) = J^{Rk) , (75) 

whose numerical solution gives the value of kR for a given s and m, and from this the 
energy levels are computed using Eg^rnj = shvpks^mj- The eigenvalues Eg^mj are defined 
by three quantum numbers: s, m, and j, where j represents the ascending order of the 
values of kR that satisfy (JTE]), for a given s and m. In the next section we give numerical 
results for the energy eigenvalues. 

5.2. The circular dot in a finite magnetic field 

Let us now see how we can adapt our formalism to address the calculation of the energy 
eigenvalues of a circular dot in a magnetic field, that is we want to study the formation 
of Landau levels in reduced geometries (amusing enough, the first calculation of Landau 
levels using the Dirac equation is as old as quantum mechanics itself [92], a result that 
was forgotten by the graphene community). Experimentally this situation has been 
realised in Ref. [93]. The cyclotron motion of bulk graphene was discussed in Ref. [94] . 
The Hamiltonian (j73ll was written for a single Dirac cone. As is shown in 



Appendix C the Hamiltonian for the two Dirac cones can be written using an additional 



quantum number k = ±, associated with the valley index, reading 
rr k ( id^ + Kdy\ 

In this section we do not use the infinite mass boundary condition, but introduce 
the zigzag type of boundary condition. This will allow us to consider the presence edge 
states [95]. As we will show in the next section, these states are always present in 
graphene quantum dots. Recalling Fig. [H one sees that at the zigzag edge only one 
type of carbon atom (either A or B) is present. The boundary condition at a zigzag 
edge with, say, only B atoms present, requires that the amplitude of the wave function 
at the A atoms to be zero, we therefore have the condition 

MR,f) = o. (77) 
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In the following we will choose the boundary condition defined by Eq. (I77l) . for which 
the wave vector is quantised as k = Zmj/R, where Zmj denotes the j-th root of the m-th 
Bessel function, Jm{zmj) = 0. 

A magnetic field B = Bcz, perpendicular to the graphene sheet, gives rise to a 

— * 

vector potential, which in polar coordinates reads A = A^e^,, and, using Gauss' theorem, 
one obtains 2'7irA^p = nr^B. Making the traditional minimal coupling of the charged 
electrons to the vector potential, the Hamiltonian has the form 



= -ihvf 











(78) 



where $o = h/e ~ 4136 T-nm^ denotes the elementary fiux quantum and — e is the 
electron charge. We now make the observation that the trial function 

renders the eigenvalue problem a one-dimensional one, with the radial Hamiltonian given 
by 



(79) 



= -ihvf 





a Km ,^ nBr 
' r <Po 







. TvBr 
' *0 



(80) 



Let us make the substitution ilf = ip^ yr (i = 1, 2) in the eigenvalue equation defined by 
the Hamiltonian (l80l) . with the radial spinor wave function having the form ip = ii'^, ^^)- 
This substitution was considered before in the exact solution of the Coulomb problem 
in the 2+1 dimensional Dirac equation [96] and also in [89]. This procedure leads to a 
more symmetric eigenproblem of the form 



ihvf 



ihvf 



r + 



Km + 1/2 TfBr 
r $0 
TcBr 



Km 



1/2 



+ K- 



(81) 



(82) 



r $0 

We want to solve the eigenproblem defined by (iSTi ) and (l82l) by diagonalising an 
Hermitian matrix. To this end let us look at the problem of a dimerised one-dimensional 



t(n-1) 



t(n-1,n) 



A B 



n-1 



n+1 



Figure 10. (colour on-line) Representation of a dimerised chain of atoms A and B. 
Within the unit cell n the hopping is t{n) and between the unit cells n and n + 1 the 
hopping is t{n, n + 1). They can be functions of the unit cell position n. 



tight-binding model, such as that represented in Fig. [TOl The relevance of this interlude 
will be apparent in a moment. The Hamiltonian for the depicted system is 

H = '^[t{n)\nA){nB\ + t{n,n + l)\nB){n + 1A\ + H. c] , (83) 

n 
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and the wave function is written as 

I*) = 5^(an|nA)+6„|nS)). (84) 

n 

The eigenvalue equation = El'if) can be reduced to the solution of the linear 

homogeneous system 

bnt{n) +bn-it{n -l,n) = anE , (85) 
ant{n) + an+it{n,n + 1) = bnE . (86) 

Introducing the simplifying notation t{n,n + 1) = t'{n), the above eigensystem reads 

bnt{n) + bn^it'in - 1) = anE , (87) 
ant{n) + an+it'{n) = bnE . (88) 

We note that Eqs. (l87ll and (l88l) pose well defined numerical problem for well behaved 
functions t{n) and t'{n). Let us now see what kind of continuous model follows from this 
lattice problem. Notice that since the model under consideration has a valence and a 
conduction bands, in case we have one electron per site the relevant energies are around 
zero. In this case, the amplitudes a„ and 6„ oscillate between positive and negative 
values within a lattice unit cell. In order to construct a well defined continuous model 
we need to subtract this oscillatory behaviour making the replacement a„ = i(— )"a„ 
and bn = {—)"'bn- This produces the set of equations 

- z[M(n) - bn-it\n - 1)] = dnE , (89) 
-i[-ant{n) + dn+it'{n)] = bnE . (90) 

Defining now T(n) = [t{n) +t'(n)]/2 and A(n) = [t{n) — t'{n)]/2 and recalling that the 
first order derivatives can be approximated by 

drd [a(r„+i) - a(r„)]/Ar , (91) 
drb^[b{rn)-b{rn-i)]/Ar. (92) 

and that Ar = R/Ni, rn = Rn/Ni a discretised position vector, with R the length of 
the chain (which will correspond latter to the radius of the dot), n = 1, . . . Ni, and A*"; 
the number of points in which the length R was discretised, we obtain 

- i[T{r)Ardrb + 2A(r)6] = dE , (93) 

- i[T{r)Ardrd - 2A(r)a] = bE . (94) 

We can thus make the following identification to the continuous model: 

nr) =1^. (95) 
±2A(r) = ^H!i±iZ? + = . (96) 



The ambiguity introduced by the ± sign in Eq. (1961) can be settled by looking at the bulk 
limit of the problem. The choice that gives the correct answer is t'{n) = T{r) + Q{r)/2 
and t{n) = T[r) — Q{r)/2. If we assume that the wall of the dot is located at n = Ni, 
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then the boundary condition (i77l) is imposed considering t/'at, = 0. In order to keep the 
problem particle-hole symmetric, we consider the case where the effective chain problem 
has Ni unit cells. 




-30 -20 -10 10 20 30 -30 -20 -10 10 20 30 -1(X) 



Figure 11. (colour on-line) Forty energy levels of a circular graphene quantum-dot 
in a finite magnetic field as function of the angular momentum quantum number m. 
The radius of the dot is i? = 70 nm, three magnetic fields were used B =1, 10, 30 T 
(from left to right), and k — 1. The two surface states (for m > 0) are represented 
using squares and circles. 

In Fig. dHwe represent the numerical solution of Eqs. ( l87l l and (l88D . It is clear that 
at small fields the bands are essentially symmetric for positive and negative m values, a 
property that comes from the fact that exchanging m by —m in the Dirac equation only 
replaces the role of ipl^ and ip^. With a finite magnetic field the situation changes. Also 
seen is the presence of a dispersive edge state. The dispersive part, occurring for positive 
m, is dependent on the Dirac point. Note the appearance of the zero energy Landau 
level upon increasing the magnetic field. The different behaviour, for large fields, shown 
by the energy levels for positive and negative m is associated to the amount of angular 
momentum induced by the magnetic field. 

Let us now discuss how to include the Coulomb interaction in the calculation. This 
is important because the screening in the dot may not be very effective and because the 
dot may be working under a regime where it has a net charge density (charged dot). The 
situation of a charged dot is represented in Fig. [121 The gate potential Vg induces either 
holes or electrons in the dot. This causes a situation where the dot is charged, since the 
neutrality case takes place when the chemical potential is at the Dirac point. If we take 
the graphene to be at a potential Vg then the charges accumulated in the metal-insulator 
interface have to be at zero potential. This means that one must use a set of images 
charges at a distance t from the interface with exactly the same spatial density of that 
formed in the graphene dot. We then have to describe the Coulomb interaction of an 
electron in the dot with both the self-consistent charge density in the graphene and its 
image underneath the metal-insulator interface |100| . The simplest way to include the 
effect of Coulomb repulsion is by using the self-consistent Hartree approximation. We 
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Figure 12. (colour on-line) Lateral view of graphene FET. When the dot is gated the 
charge accumulates at the metal-insulator (Si+-Si02) interface. The thickness of the 
insulator is t and the applied gate potential is Vg. 



analyse here the effects induced by increasing the number of electrons in the dot using 
the Hartree approximation [97j. The self-consistent Hartree potential describes, within 
a mean field approximation, the screening of charges within the dot. We assume that 
a half-filled dot is neutral, as the ionic charge compensates the electronic charge in the 
filled valence band. Away from half filling, the dot is charged. Then, an electrostatic 
potential is induced in its interior, and there is an inhomogeneous distribution of charge. 
We describe charged dots by fixing the chemical potential, and obtaining a self-consistent 
solution where all electronic states with lower energies than the Fermi energy are filled. 
From this calculation we obtain the Hartree electronic energy bands. The Hartree 
approximation should give a reasonable description when Coulomb blockade effects can 
be described as a rigid shift of the electrostatic potential within the dot [Ml IMI- Using 
the same discretisation procedure as before, the numerical equations to be solved have 
now the form 

VH{n)an + hnt{n) + &„_it'(n - 1) = an-Ej,m , (97) 
VH{n)hn + ant{n) + an+it'{n) = bnEj^^ . (98) 

where the Hartree potential in the continuum, Vnir), is given by 



r'dr'dplC{r,r',t) ^Vo—Y] / dpr'^K,{rn,r'^,t) , 



(99) 



and 



nry.t) = P^''^ - ^ P^''^ (100) 

A^r^ _l_ (j^iy — 2rr' cosy^ a/^^ + {r'Y — 2rr' cosy? + 4/!:^ 

with the parameter Vq given by Vq = (e^/47reoe) and p(r„) the electronic density at point 
Tn, computed from 

p(r„) = c/C^ K^n^j^ + bl^^jJ/rn, (101) 

jm 7^Jm,spur. 
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such that the sums over m and jm are constrained to those energy levels such that 
Em,j < Ep (note that the explicit dependence of and 6„ on m and jm has been 
introduced in Eq. llOip . where Ep is the Fermi energy measured relatively to the 
Dirac point, g is the spin and valley degeneracy, and the constant C is given by the 
normalisation condition on the disk, C = (27rAr)~^. The constraint in the j summation 
in Eq. (llOip is due to the fact that the boundary conditions introduced by the finite 
tight-binding chain fails to reproduce accurately the boundary condition ipi^r ^ 0) ^ 0, 
introducing a spurious mode, characterised by the quantum number jVrt.spur. = Ni + 1 
for m > 0; in order for sensible results to be obtained these modes have to be removed. 

We note that the integral (i99l) is well behaved since the self-interaction has been 
excluded {n ^ n'). Further the p(r„/) is independent of y?. The angular integral in Eq. 
( l99l) can be formally computed leading to 



The elliptic integral K(m) can be approximated by an analytical function |101| . which 
reduces the numerical effort. It is now clear that, due to the Hartree potential, the 
problem defined by Eqs. (I97ll and (l98l l has to be solved self-consistently. 

One should comment on the fact that for a large dot R ^ t the contribution 
from vuir) essentially vanishes and the change of the bands due to the Hartree term is 
vanishingly small. On the contrary, for small dots i? ~ t and the Hartree renormalisation 
of the electronic energy levels can be very important in the case of heavily charged dots. 

In Fig. [13] we represent the energy bands of a quantum dot of radius R = 100 nm 
on top of a silicon oxide slab of thickness t = 100 nm. We used a gate voltage of = 2 
V, which corresponds to a Fermi energy of Ep = 0.077 eV for the bulk system. The 
values of the magnetic field used were B = 1 T (top panels) and -B = 3.5 T (bottom 
panels). It is clear that the Hartree bands are renormalised by the Coulomb interaction. 
In Fig. [T4]we represent the self-consistent density, p(r), and Hartree potential, vuiji), 
for two different values of the magnetic field. The parameters are those given in the 
caption of Fig. [131 The increase in the Hartree potential upon increasing B is due to 
the increase of the number of electrons in the dot. 

We should comment that in our calculation we have not tried to keep the number 
of electrons fixed. This can easily be done, but increases the computational effort since 
the chemical potential has to be self-consistently determined. Instead, we have chosen 
to keep the Fermi energy constant, which, of course, leads to a changing in the number 
of electrons in the dot with the variation of the magnetic field. Also the population of 




(102) 



with K(m) defined as 




(103) 



Dirac electrons in graphene-based quantum wires and quantum dots 



26 




Figure 13. (colour on-line) Independent and Hartree energy levels for a spherical 
quantum dot with R = 100 nm, on top of a silicon oxide slab of t = 100 nm (e = 3.9), 
for a Fermi energy Ep — 0.077 eV (represented by a dashed line). For B = 1 T, the 
number of electrons in the dot is Ne = 166 and the magnetic length is £b = 26 nm; for 
B = 3.5 T, the number of electrons in the dot is Ne = 178 and the magnetic length is 
£b = 14 nm. The left panels are the independent energy bands; the right ones are the 
Hartree bands. The top row is for i? = 1 T. The lattice has Ni — 100. 




Figure 14. (colour on-Hne) Self-consistent electronic density (left) and Hartree 
potential VH{n) (right) for the same parameters given in the caption of Fig. [I3l 



the surface states was not included in the calculation, using a criterion of computational 
simplicity. In a future study we shall relax these two constraints. 

Another important aspect in quantum dot physics is that of confinement introduced 
by the potential creating the dot. The confinement potential can be either due to etching 
or to applied gates. In the case of dots or narrow channels described by the Schrodinger 
equation, a very popular confinement is that introduced by a parabolic potential |102| . 
since it allows a simple analytical solution. We choose a confinement potential given by 

V,onf{r)=U,{r/R)\ (104) 

which rises smoothly from the centre of the dot. The prefactor Uc is the strength of the 
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potential at the edge of the dot. As in the case of the Hartree potential, Kon/('^) enters 
in the diagonal part of the radial Hamiltonian. In Fig. [15] we give a comparison of the 
energy bands for B = 10 T. Comparing the left and the right panels of Fig. [I5] we see 
that the Landau levels become dispersive with m due to the confinement, a result also 
found for Landau levels derived from the Schrodinger equation |102| . Interestingly, we 
see that the confinement also breaks the particle-hole symmetry of the problem, a result 
found before for the ribbon problem |103| . 




-100 -50 -100 -50 -100 -50 

m m m 



Figure 15. (colour on-line) Energy spectrum of a grapliene quantum dot at -B = 10 
T. From left to right: independent particle bands, Hartree bands, and independent 
particle bands with the confinement potential l|104p . The parameters are J7c = 0.1 eV, 
Vg = 5\, Ep = 0.12 eV, and the remaining parameters are those used in Fig. [131 The 
magnetic length is £b = 8 nm and the number of electrons in the dot is TVg = 463. 



5.3. The hexagonal and circular dots at the tight-binding level 

In this subsection we want to address the question whether graphene quantum dots will 
have or not edge states, starting from the full solution of the tight-binding Hamiltonian 
([Toll . Edge states in graphene nanostructures are of particular importance since they 
can give rise to magnetism see, e.g., Ref. |104| . In order to access the low energy density 
of states of dots with physically relevant sizes (bigger than 10 nm) a Lanczos technique 
is used, since the exact diagonalisation of systems of this size becomes intractable. For 
a brief introduction to the Lanczos technique, see e.g. Ref. |105) . 

We have chosen to diagonalise dots of circular and hexagonal shape with zigzag- 
termination, such as those depicted in Fig. [HI Our numerical findings are represented 
in the Fig. [ITl It is clear that as the size of the dot grows larger the number of zero 
energy states increases, indicating the presence of zero energy-edge states. When a finite 
f is added to the Hamiltonian, the edge states become dispersive [39] and there is a 
reduction of the density of zero energy states (seen in the hexagonal dot) . 
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Figure 16. (colour on-line) Representation of a graphene crystallite on top a which 
an hexagonal or circular dots can be patterned. The size of the dot is defined by the 
length R. 




Figure 17. (colour on-Hne) Density of states at low energies for hexagonal (left) and 
circular dots (right). The insets are a zoom in of the density of states close to the 
Dirac point. Dots of several sizes are represented. 

6. Final comments 

In this work a description of the confinement of Dirac electrons in nano-wires and 
quantum dots was given. It was shown that, in principle, it is possible to localise 
electronic modes in a spatial region of a nanowire, using ap — n — p gate potential setup. 
The energy spectrum of quantum dots in a magnetic field was described taking into 
account both the effect of electron-electron interactions, at the Hartree level, and the 
effect of confining potentials. The inclusion of exchange |106l 11071 HQS) in this study can 
in principle be done. The interesting aspects about this possibility are two-fold: first, 
the exchange energy for Dirac electrons is different from that for the two-dimensional 
electron gas described by the Schrodinger equation; second, and contrary to the Hartree 
potential, the exchange energy of the full electronic system has to be considered, since 
there is not an equivalent cancellation effect to that found in the Hartree potential 
between the ion background and the valence electrons direct Coulomb energy. These 
aspects will be pursued in a follow-up study |109| . 
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Appendix A. Equations for the numerical solution of the scattering problem 

Below we give the equations that have been solve numerically when we studied the 
scattering problem by a infinite-mass wall. These are: 

• p = 0: 

2 

[(1 + izo) + ro(l + izo)] + — [(^0 + + ro{zo + i)] + 

zn 

_2 2 

Y] -. [rm{l + iZm)]+ Y] -, —— [rm{Zm + i)] ^ , 

"lodd>0 meven>0 ^ ' 

(Al) 

• p 7^ and even: 

2 

rp(l + iz^ + - — — — [(^0 + + '^ol^o + i)\ + 
lyp + IjTT 



^-^ zip — mjTT ^-^ zip + m + 1 TT 

"lodd>0 ^ ' mm,even>0 ' 

(A.2) 

with the possibility of having m = p in the meven summation 
p odd 

2 

rJ^Zp + i) + -- — \{zq + i) + ro(5o + i)\ + 

•/ , 1 1 ^ [rm(l + «5m)] + y \ ^ [rm(5m + «)] = 0, 

' 7 1 n 4- I 4- m. l-TT ' i \ — n 4- m \7r 



i(p + l + m)7r ^ i(—p + m)7r 



(A.3) 



with the possibility of having m = p in the modd summation. It is clear that we 
can truncate this set of equations to obtain a set of N equations for the coefficients 



Appendix B. General boundary conditions in a quantum dot with infinite 
mass confinement 

In this Appendix we give all the details of how to obtain the boundary condition of the 
wave function at the wall of a quantum dot, with the confinement determined by the 
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infinite mass condition. We must compute the wave function in the domain D due to a 
refiection at the boundary. We write the plane-wave inside D as 

^D=(^ j e^'^-UR j e^'f-^, (B.l) 

where from Fig. [9] we can conclude that 9i = 7r + 2a — 6*0, and ki and kf are the momenta 
of the incident and refiected waves at the boundary of the dot. In order to calculate 
ip2/i'i (which will make it possible to compute the value of B), we need to discover the 
value of the reflection coefficient, R. We can accomplish this, using the fact required 
by the Dirac equation, that the components of the spinors must be continuous at the 
boundary. 

First we solve the Dirac equation with a mass Mvp > E. For the sake of simplicity, 
we use the normal and tangential coordinates n and s given by 

n = X cos a + ysina 
s = — xsina + y cos a 

which implies that, 

dx = (dxn) dn + {dxs) 
dy = {dyu) dn + (dys) (9^ 
dx = cos adn — sin adg 
dy = sin adn + cos ads 

resulting in 

dx±tdy = {dn±tds)e^'". 
Then, for a plane wave in the domain OD (complementary to D), \E'od = 
T ^ ^ ^ gj(fc„n+fcas)^ (where T stands for the transmission coefficient). Solving the Dirac 
equation explicitly we obtain 




E~Mv% 
ihvp{q—k) 



^iks-qn ^ ^g_2) 



and chosen the solution that decays for r +00. Further we identified kn = iq and 
kg = k. Imposing the continuity of the wave functions (IB. II) and (IB. 21) at the boundary 
of the dot, one obtains 

l + R =T, 
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Replacing the value of R in Eq. f IB.ll ). the wave function in the dot reads 





which, after some simple manipulations, allows us to conclude that 

— = le 

^1 

and therefore B = 1. 

Appendix C. The Dirac equation in polar coordinates 

To treat problems with circular symmetry, the partial derivatives with respect to 
Cartesian coordinates shall be written in polar coordinates (r, ip). For the x-coordinate, 
the product rule yields dx = {dr/dx)ydr + {dip/dx)yd^, where the derivatives are taken 
for fixed y. The first derivative is obtained using r = a/x^ + y'^. The second one uses 
tanv? = y/x and thus (1/ cos^<^)9^ = —{y/x'^)dx. This gives 

dx = cos(pdr — , (CI) 



and analogously 



dy = sinv3(9r + . (C2) 



The Hamiltonian thus reads 



where we have introduced the additional quantum number k, to account for the two 
non-equivalent Dirac points. Let us now define the operators 

L+ = e'^{dr + '-d^) , (C.4) 

L_ = -e-'^'idr - '-d^) , (C.5) 

which acting on the product of a Bessel function of integer order m, Jm{kr), and a 
complex exponential, e*™''^, produce L±Jmikr)e'^"^'^ = —kJm±i{kr)e''^"^^'^'>'-^. This last 
result leads to the construction of the wave function of the free problem in the form given 
in Eq. fFTil) . In addition, the following commutators [L^,L±] = ±L± and [L+,L_] = 
allow us to interpret the operators L± as rising and lowering operators of the angular 
momentum. 

Finally we note that if we consider a ring instead of a disk it is possible to add a flux 
through the ring, introducing a vector potential = ($/27rr)e^. The full Hamiltonian 
with both a perpendicular magnetic field and the magnetic flux through the ring is given 

by 

and its numerical solution can be accommodated within the explained method. 
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